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ABSTRACT. 

In  the  case  of  nearly  incompressible  elastic  materials  the  strain  energy,  the 
shear  stress,  and  the  difference  of  normal  stresses  can  be  computed  accurately  by 
direct  methods  when  the  p-version  of  the  finite  element  method  is  used.  Compu¬ 
tation  of  the  sum  of  the  normal  stresses  requires  special  procedures.  In  this  paper 
such  procedures  are  described.  Examples  are  presented. 
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1.  INTRODUCTION. 


^  The  problem  of  computing  stresses  accurately  and  reliably  in  the  case  of  nearly 
incompressible  solids  and  in  the  incompressible  limit  of  the  compressible  Navier- 
Stokes  equations  has  received  a  great  deal  of  attention  in  recent  years.  Formu¬ 
lations  based  on  the  principle  of  virtual  work  have  been  abandoned  in  favor  of 
mixed  and  penalty  methods  in  order  to  avoid  locking  which  occurs  when  elements 
of  low  polynomial  degree  are  used.  An  overview  of  the  evolution  of  mixed  and 
penalty  methods  and  a  clear  exposition  of  the  locking  phenomenon  is  available  in 
fl|.  A  mixed  formulation  which  has  been  successfully  applied  to  elastic  and  inelas¬ 
tic  problems  is  presented  in  $\\  Mixed  formulations  have  been  used  successful iy 
in  fluid  dynamics  in  conjunction  with  elements  of  high  polynomial  degree  ,{3|. 

Unlike  the  formulations  presented  in  [1,2,3],  our  formulation  is  based  on  the 
principle  of  virtual  work.  The  solution  is  obtained  in  terms  of  displacements.  The 
sum  of  normal  stresses  is  recovered  by  post-solution  operations.  Our  approach  to 
circumventing  the  locking  problem  is  based  on  the  following: 

It  has  been  established  theoretically  and  demonstrated  numerically  that  the 
error  measured  in  energy  norm  converges  at  the  same  rate  independently  of  Pois¬ 
son’s  ratio  when  the  p- version  of  the  finite  element  method  is  used  [4].  When 
the  material  is  nearly  incompressible  (Poisson’s  ratio  is  close  to  1/2)  then  direct 
computation  of  the  sum  of  the  normal  stresses  from  the  finite  element  solution 
yields  very  inaccurate  results  because  the  volumetric  strain  is  nearly  zero  in  the 
least  squares  sense.  Small  oscillations  about  the  mean  value  occur,  due  to  errors 
of  approximation  in  the  finite  element  solution.  Direct  computation  of  the  sum 
of  the  normal  stresses  involves  multiplication  of  the  volumetric  strain  by  a  large 
number,  the  bulk  modulus,  which  magnifies  the  errors  of  approximation  [5].  For 
this  reason  the  sum  of  normal  stresses  has  to  be  computed  indirectly.  On  the  other 
hand,  the  shear  stress  and  the  difference  of  the  normal  stresses  can  be  computed 
from  the  finite  element  method  with  good  accuracy. 

The  situation  is  similar  in  fluid  dynamics  when  the  incompressible  limit  of 
the  compressible  Navier-Stokes  equations  is  approached  but  the  terminology  is 


A-,?' ' 


different:  We  speak  of  rate  of  deformation  instead  of  strain  and  pressure  instead 
of  the  sum  of  normal  stresses. 

The  indirect  method  described  herein  utilizes  the  fact  that  the  sum  of  the 
normal  stresses  satisfies  Laplace’s  equation  when  the  body  forces  are  zero  or  con¬ 
stant.  To  solve  Laplace’s  equation  it  is  necessary  to  specify  either  the  sum  of  the 
normal  stresses  or  its  derivative  with  respect  to  the  normal  along  the  boundary 
of  the  solution  domain.  Since  we  cannot  compute  these  values  from  the  finite  ele¬ 
ment  solution  directly,  we  utilize  either  specified  traction  values  or  the  equilibrium 
equation  and  the  shear  stresses  and  the  differences  of  normal  stresses  computed 
from  the  finite  element  solution. 


2.  BACKGROUND. 


We  denote  the  solution  domain  by  n  and  the  thickness  of  the  elastic  body 
by  t,  (constant),  the  displacement  vector  by  u,  and  its  Cartesian  components  by 
ux  =  uz(x, y),  uy  =  u„(x,y).  The  strain  components  corresponding  to  u  are  denoted 
by  eiu) ,  ,  'riy* .  By  definition: 

(u)  def  dtix  (u)  def  duy  (u)  d«f  dux  duy 

<9x  ’  *  ~  3y  '  dy  dx'  1  J 


In  the  case  of  plane  strain,  the  stress  components  axe  related  to  the  strain  com¬ 
ponents  by  Hooke’s  law: 


4tt)  -A  (4tt)  +  4U))  +  2G4tt) 

(2a) 

<4U)  =A  (4tt)  +  4U))  +  2G4u) 

(26) 

<4tt)  =a  (4“>  +  4“>) 

(2c) 

Uy  =<?7 

(2d) 

where 

\  Ev  c  E 

(l+i/)(l~2v)  “  2(1  +  v) 

(3) 

are  the  Lame  constants.  In  matrix  form  (2a,b,c)  are  written  as: 

{*M}=[E]{^} 

(4a) 

where 

{<r(“)}  =  {a<«>  <r<">  =  {«<,«>  e‘“>  7‘“»}T 

(46) 

and 

A  +  2  G  A  0 

[E)  =  A  A  +  2G  0  •  (4c) 

L  0  °  G 

[E\  is  a  positive  definite  matrix  and  therefore  it  can  be  written  as: 


where  [D]  is  the  diagonal  matrix  of  the  eigenvalues  of  [15]  and  [T]  is  the  matrix  of 
normalized  eigenvectors  of  [£j.  Specifically,  [D\  and  [r]  are: 


2(A  +  G)  0  O' 

‘1/V2 

l/y/2  O' 

[D\  = 

0  2G  0 

PI- 

1/V2 

-l/v/2  0 

0  0  G 

0 

0  1 

Note  that  [T\  =  [T\T  =  [Tj-1.  The  strain  energy  is: 

t-'(i)  =  \[ fa  +  rM-iff)  t.  dzdy 

The  energy  norm  of  u  is  denoted  by  ||u|[£: 


(6) 


(7a) 


Mb  =  y/m- 


(7b) 


Using  (5)  and  (6)  we  have: 

UW  =  \  ffa  [(A  +  G)(4“)  +  4U))2  +  <?(4“’  ~  4U))2  +  0[- 7<“>)2]  t,  dzdy.  (8) 


We  denote  the  exact  solution  (in  the  virtual  work  sense)  by  uEX-  To  obtain 
finite  element  approximations  to  uEX  we  construct  a  mesh  A  by  subdividing  n  into 
triangular  and/or  quadrilateral  elements,  some  or  all  of  which  may  have  curved 
sides.  The  meshes  are  so  constructed  that  any  two  elements  in  the  mesh  either 
have  a  vertex  or  an  entire  side  in  common,  or  no  points  in  common.  We  denote 
the  number  of  elements  by  Af(A).  The  functions: 

•-QPllvh  v-Qjftf,*),  (£,>?)  e  or  (£,r?)  s  n^1,  (t  =  1,2,...M(A))  (9) 

map  a  suitably  defined  standard  triangular  element  (resp.  standard  quadri¬ 
lateral  element  n^*)  into  the  »th  triangular  (resp.  quadrilateral)  element.  The  set 
of  all  functions  which  have  finite  strain  energy  on  n  and  are  polynomials  of  degree 
less  than  or  equal  p  on  fi**/,  (resp.  polynomials  of  degree  less  than  or  equal  to  p, 
supplemented  by  monomial  terms  of  degree  p  +  i  on  )  is  denoted  by  Sp(n,  A,Q) 


or  simply  5.  Any  function  in  S  is  continuous  on  fl.  We  define  the  set  of  admissible 
functions  to  be  those  functions  in  S  which  satisfy  the  kinematic  boundary  condi¬ 
tions  and  denote  them  by  S.  The  dimension  of  S  is  called  the  number  of  degrees 
of  freedom  and  is  denoted  by  N. 

The  finite  element  solution  SFE  is  that  function  from  S  which  minimizes  the 
strain  energy  of  the  error: 

U  (Hex  -  ufe)  =  min  U  (Hex  -  u)  •  (10) 

a  <5 

Denoting  «  =  uEx  -  ure,  we  have 

U(e)  =  yj[(\  +  G)  (4e>  +  4'))  3  +  (4*>  -  4*))  2  +  G  (7W)  3]  t,  dz  dy.  (11) 

Th*>  strain  energy  of  the  error  U(e)  can  be  progressively  reduced  by  constructing 
sequences  of  discretization  Slt  S2,  S3,  ...  by  successive  mesh  refinement,  increase 
of  the  polynomial  degree  of  elements,  or  both.  U(e)  approaches  zero  at  a  rate 
which  depends  on  uEx  and  the  way  in  which  the  finite  element  spaces  S2>  S3 . . . 
are  constructed.  Details  are  available  in  [6]  and  the  references  listed  therein.  Ob¬ 
taining  a  sequence  of  finite  element  solutions  corresponding  to  Si,S2, . . .  is  called 
extension.  If  the  sequence  of  finite  element  spaces  is  constructed  by  successive 
mesh  refinement  then  it  is  called  h-extension,  if  the  mesh  is  fixed  and  the  polyno¬ 
mial  degree  is  increased  then  it  is  called  p-extension  and  if  the  two  are  combined 
then  it  is  called  hp-extension.  In  extensions  the  number  of  degrees  of  freedom  is 
progressively  increased. 

It  has  been  shown  theoretically  and  demonstrated  numerically  that  in  p- 
extensions  U(e)  0  at  a  rate  which  is  independent  of  u,  even  if  1/  is  very  close 
to  its  limiting  value  of  1/2  [4].  In  the  case  of  h-extensions  and  small  p  the  rate 
of  convergence  slows  very  substantially  as  v  -*  1/2.  The  reason  for  this  can  be 
seen  from  (11):  As  v  1/2,  X  -*  00  and  the  first  term  of  the  integrand  in  (11)  can 
be  viewed  as  a  penalty  term,  i.e.  a  constraint  which  requires  that  (e*  1  +  4"’)  be 

very  nearly  zero  in  the  least  squares  sense.  One  constraint  is  associated  with  each 
element,  thereby  reducing  the  number  of  degrees  of  freedom  by  one  per  element. 
When  h-extension  is  used  and  p  is  low  then  the  number  of  constraints  increases 
at  almost  the  same  rate  as  N  and  the  rate  of  convergence  slows.  In  the  case  of 
p-extension  the  number  of  constraints  does  not  grow,  hence  nearly  incompressible 
materials  can  be  analyzed  without  significant  loss  in  effectiveness. 


3.  COMPUTATION  OF  STRESSES. 


We  now  define  the  root-mean-square  stress  measure  S( u)  as  follows: 

m  -  \jl  ffn  ’ +(»!*') 2 + ('];’)  i «.  *  *  (i2) 


where  V  is  the  volume.  Using  (4a,b.c)  and  (5)  S2(u)  can  be  written  in  terms  of  the 
strains  as  follows: 

S2(u)  =  i  ||n[2(A  +  G)2(4“>+£<“))3+2G2(e(“>-e«“>)2  +  G2(7i“»)2]  t.dxdy.  (13) 

On  comparing  S2(d)  with  l/(u)  in  (8)  we  see  that  the  two  expressions  are  similar. 
In  fact,  we  can  write: 


•$?(«)  d=  Jfn2(X  +  G ^  (4*e)  +  4*’)  tz  dxdy  =  W  Jfn  (a*e)  +  4*')  tx  dx  dy 


(14a) 


<-(A  +  G)U& 


that  is,  the  square  of  the  error  in  the  sum  of  the  normal  stresses  integrated  over  the 
volume  is  bounded  by  the  strain  energy  of  the  error  U(e)  multiplied  by  a  constant 
which  depends  on  Poisson’s  'atio  and  goes  to  infinity  as  u  -*  1/2.  On  the  other 
hand,  the  error  in  the  differences  of  normal  stresses  and  shear  stress  can  be  made 
arbitrarily  small  even  if  v  — ■  1/2: 


Sl{i i)  =  £  ffo  2 G3  -  4e>)3  t.  dxdy  = -L  (<4«>  -  <r<'))2  t.  dxdy<±  GU(e)  (14b) 

and 

SI  W  d=  vJJqG2  (*')’  t.d*dy=±J ^  (r<'»)2  1,  dxdy  <  | GU(i ).  (14c) 

This  analysis  indicates  that  we  can  expect  good  approximations  (in  the  mean)  to 
uz  -  av  and  rxv  from  the  finite  element  solution  via  the  constitutive  relationships 
(2a-d)  but  not  to  a x+oy.  Therefore  we  must  avoid  using  the  constitutive  equations 
in  the  computation  of  the  sum  of  normal  stresses.  To  do  this  we  utilize  the  fact 


that  when  the  body  forces  are  zero  or  constant  then  the  sum  of  normal  stresses 
satisfies  Laplace’s  equation: 

A{£T,  +  (Ty)  =  0.  (15) 


For  details  see,  for  example,  [7,8].  Obvious  modifications  can  be  made  for  the 
case  of  general  body  forces.  The  boundary  conditions  are  established  from  the 
equilibrium  equations  using  the  fact  ai and  rlXyrz)  can  be  computed  from 
the  finite  element  solution  with  sufficient  accuracy.  We  denote  typical  boundary 
segments  of  the  solution  domain  fl  by  T(») ,  (i  =  1,2, 3,4)  and  we  denote  the  positive 
(outward)  unit  normal  by  n,  the  positive  unit  tangent  by  t.  Four  cases  are  possible: 

(1)  *n,  rnt  are  given  on  I*i.  Since  the  sum  of  normal  stresses  is  invariant  with 
respect  to  coordinate  transformation: 


<TX  +  <Jy  =  (T„  +  <Tt  ~  ~  (fn'*  1  “  +  2crj,ri 


where  the  superscript  (rj)  indicates  that  the  function,  in  this  case  an ,  is  known 
from  the  specified  boundary  condition. 

(2)  u,  are  given  on  r2.  The  procedure  is  the  same  a s  in  case  (1)  with  r2  replacing 
Ti  in  (16).  An  important  special  case  is  the  condition  of  antisymmetry,  i.  e. 
<rn  =  0,  14  =  0.  In  this  case  crx  +  ay  =  0. 

(3)  «n,  rnt  are  given  on  r3.  In  this  case  the  normal  derivative  of  <jx+<ry  is  computed 


from: 


dWx+<ry)  =  d(cn  +  <Jt )  =  3(al“,,',)  -a{tu,‘))  _ 

dn  dn  dn  dt  3 


where  we  used  the  equilibrium  equation: 


9cn  dTnt_  _ 
dn  dt 


When  boundary  conditions  of  type  (1)  or  (2)  are  imposed  on  part  of  the 
boundary  then  C3  =  0.  When  boundary  conditions  of  type  (3)  are  prescribed 
on  the  entire  boundary  then  C3  is  a  constant  determined  so  that  equilibrium 
is  satisfied  by  the  boundary  conditions: 


-1- 


3 


An  important  special  case  is  the  condition  of  symmetry,  i.  e.  un  =  0,  rnt  =  0.  In 
this  case  d(<rx  +  ay)ldn  =  0. 

(4)  «n)  ut  are  given  on  r4: 

dfa  +  cry)  a(an+at)  _ 

dn  dn  dn  dt  4- 

When  boundary  conditions  of  type  (1)  or  (2)  are  imposed  on  part  of  the 
boundary  then  C4  =  0.  When  boundary  conditions  of  type  (4)  are  prescribed 
on  the  entire  boundary  then  C4  is  a  constant  determined  so  that  equilibrium 
is  satisfied  by  the  boundary  conditions: 


-/( 


d(o\rm'  -gj* 

dn 


da  +  C4  £  da  =  0. 


When  only  boundary  conditions  of  type  (3)  and  (4)  are  imposed  on  the  entire 
boundary  then  C3  =  C4  =  C  and  C  is  determined  analogously  to  (19)  and  (21). 

Thus  in  cases  (1)  and  (2)  the  sum  of  the  normal  stresses  is  specified  on  the 
boundary  (Dirichlet  condition),  in  causes  (3)  and  (4)  the  normal  derivative  of  the 
sum  of  the  principal  stresses  is  specified  (Neumann  condition).  Either  the  given 
boundary  condition  or  the  equilibrium  equation  is  used  to  avoid  computation  of 
(cr,+cry)  or  its  normal  derivative  from  the  finite  element  solution  via  the  constitutive 
equations.  When  only  cases  (3)  and  (4)  are  specified  along  the  entire  boundary 
then  the  pressure  can  be  determined  only  up  to  an  arbitrary  constant. 

The  stresses  computed  from  the  finite  element  solution  are  not  continuous  but 
the  sum  of  the  normal  stresses  specified  on  the  boundary  in  cases  (1)  and  (2)  has 
to  be  continuous.  For  this  reason  some  averaging  of  the  values  at  the  nodes  is 
necessary. 

3.1  Points  of  stress  singularity. 

When  one  or  more  points  on  the  boundary  are  points  of  stress  singularity  then 
the  procedure  just  described  has  to  be  modified  as  follows:  In  the  neighborhood 
of  singular  points  the  exact  solution  can  be  written  in  the  form: 


♦<(*)•  r<r0 


where  r,  9  are  polar  coordinates  centered  on  the  singular  point;  the  coefficients  Ai 
are  called  generalized  stress  intensity  factors;  A,  and  V>,(0)  are  determined  from  the 
conditions  that  uEx  must  satisfy  the  equations  of  equilibrium  and  the  boundary 
conditions  in  the  neighborhood  of  the  singular  points;  r0  >  0  is  the  radius  of 
convergence  of  the  infinite  series  (22). 

The  coefficients  Ai  can  be  computed  from  the  finite  element  solution  by  a 
method  so  that  the  accuracy  of  Ai  depends  only  on  the  accuracy  of  the  finite 
element  solution  HFe  measured  in  energy  norm  and  the  accuracy  of  the  finite 
element  solution  to  an  auxiliary  problem  wfe ,  also  measured  in  energy  norm, 
which  differs  from  the  original  problem  only  in  loading.  Thus: 

<  ||u£X  -  -  ti'FsIU  (23) 

The  function  wfe  is  not  actually  computed,  it  serves  theoretical  purposes  only.  In 
general  the  errors  ||uex  -  Hfe\\e  and  ||u7bx  -  wj-eIIb  are  of  comparable  magnitude, 
hence  the  error  in  the  computed  values  of  Ai  is  comparable  with  the  error  in 
strain  energy.  For  details  see  references  {9-12].  Since  the  accuracy  of  Ai  depends 
only  on  the  accuracy  of  solutions  measured  in  energy  norm  which,  in  the  case 
of  p-extensions,  is  not  sensitive  to  Poisson’s  ratio,  it  is  possible  to  compute  Ai 
accurately  independently  of  Poisson’s  ratio. 

Once  Ai  are  known,  the  loading  corresponding  to  the  stress  singular  term  can 
be  subtracted  from  the  applied  loading  and  the  problem  solved  as  a  smooth  prob¬ 
lem.  Alternatively  a  small  neighborhood  of  the  singular  point  can  be  “removed” 
from  the  domain  and  either  the  displacements  or  tractions  corresponding  to  (22) 
imposed  on  the  boundary  which  defines  the  small  neighborhood.  An  example  is 
given  in  Section  5. 


4.  EXAMPLE:  RIGID  CIRCULAR  INCLUSION  IN  AN  ELASTIC  PLATE. 

We  will  investigate  the  case  of  a  rigid  circular  inclusion  in  an  infinite  plate, 
subjected  to  unidirectional  tension.  Plane  strain  conditions  are  assumed.  The 
notation  is  shown  in  Fig.  1.  This  example  is  typical  of  cases  where  all  derivatives 
of  the  exact  solution  are  continuous  and  bounded  on  the  entire  solution  domain, 
including  the  boundary  of  the  solution  domain,  but  one  or  more  singular  points 
lie  outside  of  the  solution  domain. 


Fig.  1.  Notation. 


4.1.  The  exact  solution. 


In  this  case  the  exact  solution  is  available  [8].  The  exact  displacement  compo¬ 


nents  are: 


Ur  =  —  1)  r2  +  2 7 a2  +  1)  a2  +  2  r2  +  -  ^ -  cos2flj 


u»  =  -  ~Z  \&  (k  -  1)  a2  +  2  r2  - 


sin  2  8 


and  the  exact  stress  components  are: 


& oo  /,  £a2  35a4  \  . 

r-«  =  "  ~  (l  +  7T  +  — )  am2B 


TOW 


where  k,  p,  7,  6  axe  constants  which  depend  on  Poisson’s  ratio  u  only.  In  the  case 
of  plane  strain : 


2  2-4  v  '  1 

*  =  3-4,;  1  =  — — ; 


The  center  of  the  rigid  circular  inclusion  is  a  singular  point. 


4.2.  Finite  element  solutions. 

Finite  element  solutions  were  obtained  by  means  of  the  computer  program 
PROBE  [13].  The  solution  domain  is  defined  by  a  =  1,  w  -  b  »  4  in  Fig.  1.  The 
domain  and  two  finite  element  meshes,  a  four-element  mesh  and  an  eight-element 
mesh  are  shown  in  Fig.  2. 


D  CD  C 


Fig.  2.  Solution  domain  and  finite  element  meshes. 

The  boundary  conditions  are  as  follows:  Along  the  circular  arc  AE  both  dis¬ 
placement  components  are  zero.  Along  the  symmetry  lines  AB  and  DE  the  normal 
displacement  component  and  the  shear  stress  are  zero.  Along  boundaries  BC  and 
CD  the  tractions  computed  from  (25.a,b,c)  are  imposed.  The  corresponding  load 
vectors  were  computed  by  numerical  quadrature  using  twelve  Gauss  points  per 
element  side. 

4.2.1.  Four-element  mesh,  0  <  v  <  0.4999999. 

The  exact  value  of  the  strain  energy  for  v  =  0.4999,  computed  from  (24a, b) 
and  (25a, b,c),  is  5.79590408  a^a^t^/E.  The  number  of  degrees  of  freedom  N;  the 


strain  energy  values  computed  from  the  finite  element  solutions  corresponding 
to  p  ranging  from  1  to  8;  the  estimated  algebraic  rate  of  convergence  2/3,  and 
the  estimated  and  true  relative  errors  in  energy  norm  are  given  in  Table  1.  The 
estimated  algebraic  rate  of  convergence  20  is  an  estimate  of  the  absolute  value  of 
the  slope  of  the  log  U(e)  vs.  log  JV  curve.  The  estimated  relative  errors  in  energy 
norm,  defined  by 


Me 


def  . nn  jlggX  -  «rsl|g 

WZexWb 


(27) 


were  computed  by  the  method  described  in  [13,14].  We  see  that  the  relative  error 
in  energy  norm  is  under  one  percent  at  p  =  8. 


Table  1.  Convergence  of  the  strain  energy. 
Estimated  and  true  relative  errors  in  energy  norm. 
Four-element  mesh,  v  =  0.4999. 


p 

N 

U(ufe)E 

Est.'d 

2/3 

Est.'d 

Me  (%) 

True 

Me  {%) 

1 

8 

0.35682 

— 

96.87 

96.87 

2 

24 

3.15500 

0.66 

67.50 

67.50 

3 

40 

5.16100 

2.79 

33.10 

33.10 

4 

64 

5.66227 

3.31 

15.19 

15.18 

5 

96 

5.75195 

2.74 

8.72 

8.71 

6 

136 

5.78552 

4.12 

4.26 

4.23 

7 

184 

5.79388 

5.27 

1.92 

1.87 

8 

240 

5.79550 

5.27 

0.95 

0.84 

OO 

OO 

5.79590 

OO 

0 

0 

The  relationship  between  relative  error  in  energy  norm  and  the  polynomial 
degree  of  approximation  for  other  values  of  Poisson’s  ratio  are  shown  in  Fig.  3. 
Note  that  the  relative  error  in  energy  norm  is  plotted  on  a  logarithmic  scale  against 
the  square  root  of  the  number  of  degrees  of  freedom  N.  The  reason  for  this  is  that 
when  the  exact  solution  is  analytic  on  the  solution  domain,  including  its  boundary, 
then  the  rate  of  convergence  of  p-extensions  is  exponential: 


Me  < 


k 

ex  pilVN) 


(28) 


where  k  and  7  are  positive  constants,  [6].  Constant  k  depends  on  Poisson’s  ratio 
but  7  does  not.  This  is  demonstrated  by  the  numerical  results  summarized  in  Fig. 


3  which  show  that  the  relative  error  in  energy  norm  depends  on  u  but  its  rate  of 
change  with  respect  to  \fN  is  constant. 


0.4999999 

0.499999 

0.49999 

0.4999 


Fig.  3.  True  relative  error  in  energy  norm  vs.  y/N  for  various  Poisson’s  ratios. 

Four-element  mesh. 


4.2.2.  Eight-element  mesh,  v  =  0.4999999. 

Finite  element  solutions  were  also  obtained  with  the  eight-element  mesh  shown 
in  Fig.  2b  for  u  =  0.4999999.  The  exact  value  of  the  strain  energy  in  this  case  is 
is  5.79517057  a^afti/E.  Convergence  of  the  strain  energy  for  p  ranging  from  1  to  8 
and  the  estimated  and  true  relative  errors  in  energy  norm  are  shown  in  Table  2. 
Once  again  the  relative  error  in  energy  norm  is  under  one  percent  at  p  =  8. 

4.3.  Computation  of  the  sum  of  the  normal  stresses. 

The  sum  of  the  normal  stresses  was  computed  with  the  four-element  mesh 


Table  2.  Convergence  of  the  strain  energy. 
Estimated  and  true  relative  errors  in  energy  norm. 
Eight-element  mesh,  v  =  0.4999999. 


p 

N 

U(ufe)E 

Est.'d 

2$ 

Est.'d 
Me  (%) 

Thie 

Ms  (%) 

1 

16 

0.00029 

0.00 

100.00 

100.00 

2 

48 

0.03364 

0.01 

99.71 

99.71 

3 

80 

2.69231 

1.21 

73.17 

73.17 

4 

128 

5.52142 

5.17 

21.74 

21.73 

5 

192 

5.71015 

2.88 

12.12 

12.11 

6 

272 

5.78156 

5.25 

4.85 

4.85 

7 

368 

5.79393 

7.79 

1.49 

1.46 

8 

480 

5.79506 

7.79 

0.53 

0.44 

oo 

OO 

5.79517 

OO 

0 

0 

shown  in  Fig.  2a  for  v  =  0.4999  and  with  the  eight-element  mesh  shown  in  Fig.  2b 
for  v  =  0.4999999. 

Along  the  circular  segment  AE  the  normal  derivative  of  the  sum  of  the  normal 
stresses  was  computed  from  (20).  Along  the  symmetry  lines  AB,  DE  the  natural 
boundary  condition  d{az  +  ay)/dn  =  o  was  imposed  and  along  boundaries  BC,  CD 
<Tx+ay  was  computed  from  (16).  In  order  to  utilize  existing  capabilities  of  the  com¬ 
puter  program  PROBE  as  much  as  possible,  the  normal  derivative  was  computed 
numerically  at  nine  equally  spaced  points  along  each  of  the  two  curved  element 
boundaries.  The  normal  derivatives  were  then  determined  at  12  Gauss  points  by 
interpolation.  The  load  vector  terms  were  computed  by  Gaussian  quadrature,  us¬ 
ing  12  quadrature  points.  With  appropriate  modifications  in  the  program,  transfer 
of  boundary  conditions  from  the  solution  of  the  elasticity  problem  to  the  Laplace 
problem  can  be  automated  so  that  no  loss  in  accuracy  is  incurred  in  the  transfer 
itself. 

The  essential  boundary  conditions  were  handled  by  computing  ( ax  +  a,,)  at  ten 
equally  spaced  points  along  each  element  boundary  on  AB  and  BC,  including  the 
vertices.  At  the  vertices  the  computed  values  were  averaged.  For  example,  in  the 
case  of  the  four-element  mesh  (u  =  0.4999)  and  p  =  8,  at  vertex  C: 

(<r,  +  av)c  =  (0.9980777  +  1.0019223)  =  1.0000000  (Too  (29) 

it 

which  agrees  to  at  least  eight  significant  digits  with  the  exact  value  of  (cr*  +  <ry). 
Least  squares  approximation  was  used  to  obtain  coefficients  of  the  basis  functions 


A  EXACT 
o  EXTRACTION 
a  DIRECT  COMPUTATION 


.0  10.0  20.0  30.0  UO.O  50.0  60.0  70.0  80.0  90.0 

ANGLE  (DEGREES) 


Fig.  4a.  The  sum  of  normal  stresses  along  the  rigid  inclusion. 
Four-element  mesh,  p  =  8,  u  =  0.4999. 
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Fig.  4b.  The  sum  of  normal  stresses  along  the  rigid  inclusion. 
Eight-element  mesh,  p  =  8,  u  =  0.4999999. 


along  the  element  boundaries  with  the  restriction  that  the  values  at  the  vertices 
were  fitted  exactly. 

The  sum  of  the  normal  stresses  along  the  rigid  inclusion,  computed  directly 
from  the  finite  element  solution,  and  by  the  extraction  procedure  described  in  this 
paper,  are  plotted  along  with  their  exact  values  in  Fig.  4a  and  Fig.  4b.  The  poly¬ 
nomial  degree  8  was  used  in  solving  both  the  elasticity  problem  and  the  Laplace 
problem.  The  estimated  and  true  relative  errors  in  energy  norm  for  the  elasticity 
problem  are  in  Tables  1  and  2.  In  the  case  of  the  Laplace  problem  the  estimated 
relative  errors  in  energy  norm  were  as  follows:  Four-element  mesh,  p  =  8,  N  =  136, 
(er)s  =  0.10%.  Eight-element  mesh,  p  =  8,  N  =  2r%  ( er)E  =  0.04%. 

The  results  shown  in  Figures  4a, b  are  indicative  of  the  quality  of  approximation 
obtained  for  the  entire  domain.  A  contour  plot  for  the  sum  of  normal  stresses, 
computed  with  the  eight-element  mesh,  v  -  0.4999999,  is  shown  in  Fig.  5. 


Fig.  5.  Contour  plot  of  the  sum  of  normal  stresses. 

Eight-element  mesh,  p  =  8,  v  =  0.4999999.  Contour  interval:  0.25 
B  :  -0.75 a„,E:  0,1:  1.0O<t„,  M  :  2.OO0-OO,  P  :  2.75er0O. 

Equations  (14b, c)  show  that  ax-cy  and  rxy  converge  in  the  least  squares  sense  at 
the  same  rate  as  the  error  in  energy  norm  does,  which  is  exponential.  Exponential 
convergence  can  be  proven  for  pointwise  values  of  ax  -  ay,  rxy  and  their  derivatives 
also,  but  the  proof  is  beyond  the  scope  of  this  paper.  This  example  demonstrates 
that  control  of  errors  in  boundary  data  does  not  require  large  computational  effort. 


5.  EXAMPLE:  L-SHAPED  DOMAIN. 

This  example  is  typical  of  cases  where,  due  to  reentrant  corners,  sudden 
changes  in  boundary  conditions,  material  properties  or  loading,  one  or  more  sin¬ 
gular  points  lie  on  the  boundary  of  the  solution  domain. 


C  B 


Fig.  6.  L-shaped  domain.  Notation. 


5.1.  The  exact  solution. 

The  L-shaped  plane  elastic  body,  shown  in  Fig.  6,  is  loaded  by  tractions  which 
correspond  to  the  following  exact  displacement  field: 

ux=4  r*‘  [(*  -  <?i(Ai  +  !))  cos  w  -  A1  cos(Ai  -  2W  (30a) 

2  G 

^  rAl  [(«  +  <3X(AX  +  1))  sin \l8  +  Ax  sin^  -  2)*|  (306) 

2  G 

where  Ai  is  a  constant,  analogous  to  the  mode  1  stress  intensity  factor  in  linear 
elastic  fracture  mechanics;  Ai  =  0.544483737,  Qi  ~  0.543075579  are  constants  deter¬ 
mined  so  that  the  solution  satisfies  the  equations  of  equilibrium  and  the  stress 
free  boundary  conditions  on  the  reentrant  edges.  G  =  E/ 2(1  +  v)  is  the  modulus  of 
rigidity;  k  is  a  constant  which  depends  on  Poisson  s  ratio  and  whether  plane  stress 
or  plane  strain  conditions  are  assumed.  In  this  example  plane  strain  conditions 
are  assumed  hence  k  =  3  —  4 v.  The  stress  components  are: 

-17- 


(31a) 
(31  b) 
(31c) 


a,  =  Ax  A!  rAi_1  ((2  -  g^Ax  +  l))  cosfAx  -  1)  9  -  (Ax  -  1)  cosfAx  -  3)  9\ 

<Ty  =  Ai  Ai  rXl~1  [(2  +  gi(Ai  +  l))  cos(Ax  —  1)  9  +  (Ai  —  1)  cos(Aj  -  3)  9) 
rxy  =  Ai  Ai  rXl_1  ((Ax  -  1)  sin(Ai  -  3)  9  +  gi(Ax  +  1)  sin(Aj  -  1)  9} . 

The  boundary  conditions  are  as  follows:  Along  the  reentrant  edges  AB  and  FA  the 
normal  and  shear  stresses  are  zero;  along  edges  BC,  CD,  DE  and  EF  the  normal 
and  shear  stresses  corresponding  to  (28a, b,c)  are  prescribed.  These  tractions  sat¬ 
isfy  equilibrium  hence  only  rigid  body  constraints  have  to  be  applied.  The  exact 
value  of  the  strain  energy  is: 


U(uBx)  =  2.35967  Afa^HjE.  (32) 

Additional  information  concerning  this  model  problem  is  available  in  [12,14]. 

5.2.  Finite  element  solutions. 

In  this  example  the  eighteen-element  mesh  shown  in  Fig.  7  and  u  =  0.4999999  was 
used.  The  coefficient  Ai  was  determined  by  the  cutoff  function  method  described 
in  [12,13].  The  results  of  finite  element  solutions  corresponding  to  p  ranging  from 
1  to  8  are  given  in  Table  3. 


Table  3.  Convergence  of  the  strain  energy  and  the  computed  value 
of  the  generalized  stress  intensity  factor  Ax. 
Eighteen-element  mesh,  v  =  0.4999999. 


p 

N 

U($fe) 

U(ubx) 

(Ax)fb 

(^i)sx 

1 

41 

0.00184 

3.06910 

2 

119 

0.42400 

25.27936 

3 

209 

0.84901 

2.48314 

4 

335 

0.94668 

1.30447 

5 

497 

0.99538 

1.02696 

6 

695 

0.99761 

1.01301 

7 

929 

0.99949 

1.00296 

8 

1199 

0.99974 

1.00105 

These  results  show  that  the  errors  in  strain  energy  and  Ax  are  of  comparable 
magnitude  and  Ax  converges  at  approximately  the  same  rate  as  the  strain  energy 
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Fig.  8.  L-shaped  domain.  Sum  of  normal  stresses  at  r  =  o.la  (p 
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does.  This  is  consistent  with  the  theoretical  estimate  for  extraction  procedures 

(23).  In  general,  data  computed  by  extraction  procedures  exhibit  superconver¬ 
gence. 

5.3.  Computation  of  the  sum  of  normal  stresses. 

In  the  close  neighborhood  of  the  singular  point  the  exact  solution  is  approx¬ 
imated  by  the  expression  (22)  with  computed  by  the  cutoff  function  method 
described  in  [12].  Thus  the  error  in  the  displacements  and  stresses,  and  any  other 
quantity  computed  from  (22),  is  controlled  by  the  error  in  the  coefficients  A,.  In 
this  model  problem  the  exact  solution  is  the  first  term  in  (22)  hence  we  need  not 
be  concerned  with  the  question  of  where  to  truncate  the  infinite  series  (22).  We 
note,  however,  that  the  coefficient  of  the  second  term,  corresponding  to  the  first 
antisymmetric  mode  of  deformation  (A2  =  0.90853),  computed  from  the  finite  ele¬ 
ment  solution  of  p  =  8  by  the  cutoff  function  method,  is  2.38  x  io~8  From 

Table  3  we  see  that  At  is  accurate  to  0.1  percent,  hence  all  stress  components  are 
also  accurate  to  0.1  percent. 

In  order  to  simulate  the  case  where  the  stresses  are  computed  from  (22)  only 
in  the  immediate  neighborhood  of  the  singular  point  but  elsewhere  the  method 
described  in  Section  4  is  used,  we  removed  from  the  solution  domain  the  circular 
sector  of  radius  0.0225  a,  defined  by  the  inner  circle  in  Fig.  7.  We  computed  the 
sum  of  normal  stresses  corresponding  to  (22)  along  the  boundary  of  this  sector. 
We  imposed  the  symmetry  boundary  condition  along  the  x-axis  (see  Fig.  6)  and 
along  the  other  boundary  segments  we  computed  the  sum  of  normal  stresses  from 
(16).  The  computed  sum  of  normal  stresses  along  a  circle  of  radius  o.la  is  shown 
in  Fig.  8.  We  see  that  the  accuracy  is  about  the  same  as  in  the  case  of  smooth 
solutions  discussed  in  Section  4. 

The  results  of  this  example  are  typical  of  the  performance  of  hp-extensions 
where  the  error  in  energy  norm  is  of  the  order:  exp(— ■ y  TV1/3) ,  7  >  0.  Once  again, 
it  can  be  proven  that  convergence  of  crx  -  <r„,  rIW  and  their  derivatives,  computed 
at  points  chosen  independently  of  the  mesh  refinement,  has  the  same  exponential 
character.  This  explains  the  high  accuracy  of  the  results  in  this  example. 
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6.  NOTE  ON  THE  INCOMPRESSIBLE  LIMIT  OF  THE 
COMPRESSIBLE  NAVIER-STOKES  EQUATIONS. 


The  mathematical  formulation  of  compressible  Navier-Stohes  equations  is  very 
different  from  the  mathematical  formulation  of  the  linear  elasticity  problem.  Nev¬ 
ertheless,  in  the  case  of  viscous  fluids  which  satisfy  certain  restrictions  on  the  first 
and  second  viscosity  coefficients,  the  limiting  case  with  respect  to  the  Mach  num¬ 
ber  approaching  zero  is  exactly  analogous  to  the  limiting  case  of  incompressible 
elasticity  with  respect  to  Poisson’s  ratio  approaching  1/2.  Precise  statement  of 
the  necessary  conditio;"*  and  rigorous  proof  is  available  in  [15]. 
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7.  SUMMARY  AND  CONCLUSIONS. 

We  have  described  a  method  for  the  computation  of  stress  components  in 
the  case  of  neaxly  incompressible  isotropic  elastic  materials.  We  have  shown  that 
in  those  problems  where  the  exact  solution  is  smooth  the  shear  stresses  and  the 
differences  of  normal  stresses  can  be  computed  directly  from  the  finite  element 
solution,  provided  that  the  error  measured  in  energy  norm  is  sufficiently  small.  For 
this  class  of  problems  p-extensions  provide  efficient  means  for  reducing  the  errors 
of  approximation  and  the  performance  of  p-extensions  is  not  sensitive  to  Poisson’s 
ratio.  Thus  yield  criteria  which  are  independent  of  the  sum  of  normal  stresses, 
such  as  the  maximum  shear  stress  criterion,  can  be  applied  directly,  independently 
of  Poisson’s  ratio. 

When  Poisson’s  ratio  is  close  to  1/2  then  the  sum  of  the  normal  stresses  has 
to  be  determined  indirectly.  When  the  body  forces  are  zero  or  constant  then 
this  involves  solution  of  the  Laplace  equation,  the  boundary  conditions  of  which 
are  determined  from  the  finite  element  solution  of  the  elasticity  problem  and  the 
specified  boundary  tractions. 

In  those  problems  where  one  or  more  points  of  stress  singularity  occur  on  the 
boundary  of  the  solution  domain  the  procedure  is  similar  to  the  case  of  smooth 
solution  but  in  the  neighborhood  of  singular  points  the  approximation  is  in  terms 
of  the  coefficients  of  an  asymptotic  expansion  which  can  be  computed  by  methods 
described  in  [7-10].  Important  advantages  of  these  methods  are  that  they  are  very 
efficient  and  their  accuracy  is  not  sensitive  to  Poisson’s  ratio. 

In  the  method  described  herein  the  solution  for  the  sum  of  normal  stresses  is 
decoupled  from  the  solution  of  the  displacement  components  in  the  sense  that  first 
the  displacement  components  are  approximated,  independently  of  the  sum  of  the 
normal  stresses,  then  the  sum  of  normal  stresses  is  approximated  in  a  separate 
solution  step.  In  mixed  formulations  such  decoupling  generally  does  not  occur. 
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